home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_2 / muller.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  1.4 KB  |  62 lines  |  [MATF/MATL]
  1. function [p2,y2,err,P] = muller(f,p0,p1,p2,delta,epsilon,max1)
  2. % [p2,y2,err] = muller(f,p0,p1,p2,delta,epsilon,max1)
  3. % [p2,y2,err,P] = muller(f,p0,p1,p2,delta,epsilon,max1)
  4. % Muller's method is used to locate a root.
  5. % f is the function, input.
  6. % p0 is a starting point, input.
  7. % p1 is a starting point, input.
  8. % p2 is a starting point, input.
  9. % delta is the tolerance for p2, input.
  10. % epsilon is the tolerance for y2, input.
  11. % max1 is the maximum number of iterations, input.
  12. % p2 is the root, output.
  13. % y2 is the function value, output.
  14. % err is the error estimate for p2, output.
  15. % P is the is the vector of iterations, output.
  16. P(1) = p0;
  17. P(2) = p1;
  18. P(3) = p2;
  19. y0  = feval(f,p0);
  20. y1  = feval(f,p1);
  21. y2  = feval(f,p2);
  22. for k=1:max1,
  23.   h0 = p0 - p2;
  24.   h1 = p1 - p2;
  25.   c  = y2;
  26.   e0 = y0 - c;
  27.   e1 = y1 - c;
  28.   det1 = h0*h1*(h0-h1);
  29.   a  = (e0*h1 - h0*e1)/det1;
  30.   b  = (h0^2*e1 - h1^2*e0)/det1;
  31.   if  b^2 > 4*a*c,
  32.     disc = sqrt(b^2 - 4*a*c);
  33.   else
  34.     disc = 0;
  35.   end
  36.   if b < 0, disc = - disc; end
  37.   z = - 2*c/(b + disc);
  38.   p3 = p2 + z;
  39.   if  abs(p3-p1) < abs(p3-p0),
  40.     u = p1;
  41.     p1 = p0;
  42.     p0 = u;
  43.     v = y1;
  44.     y1 = y0;
  45.     y0 = v;
  46.   end
  47.   if  abs(p3-p2) < abs(p3-p1),
  48.     u = p2;
  49.     p2 = p1;
  50.     p1 = u;
  51.     v = y2;
  52.     y2 = y1;
  53.     y1 = v;   
  54.   end
  55.   p2 = p3;
  56.   y2 = feval(f,p2);
  57.   P = [P,p2];
  58.   err = abs(z);
  59.   relerr = err/(abs(p3)+eps);
  60.   if (err<delta)|(relerr<delta)|(abs(y1)<epsilon), break, end
  61. end
  62.